********************************************************************************************
**** TREATTEST
**** This program compares means and equality of means
********************************************************************************************

	capt prog drop treattest

	program treattest, eclass

	syntax varlist [if] [in], [ * ]
	*marksample touse
	*markout `touse'
	tempname m_mu_b m_mu_1_b m_mu_2_b m_sd_b m_sd_1_b m_sd_2_b m_coef_c_b m_d_se_c_b m_d_p_c_b cm_coef_c_b cm_d_se_c_b cm_d_p_c_b m_mu_f m_mu_1_f m_mu_2_f m_sd_f m_sd_1_f m_sd_2_f m_coef_c_f m_d_se_c_f m_d_p_c_f ivm_coef_c_f ivm_d_se_c_f ivm_d_p_c_f cm_coef_c_f cm_d_se_c_f cm_d_p_c_f  obs                       
                           
	foreach var of local varlist {
		
		* Treatment 
				
		qui ttest `var' if round == 2010, by(mother)
		
			mat `m_mu_b' = nullmat(`m_mu_b'), (r(mu_1)*r(N_1)+r(mu_2)*r(N_2))/(r(N_1)+r(N_2))
			mat `m_mu_1_b' = nullmat(`m_mu_1_b'), r(mu_1)
			mat `m_mu_2_b' = nullmat(`m_mu_2_b'), r(mu_2)
			mat `m_sd_b' = nullmat(`m_sd_b'), r(sd)
			mat `m_sd_1_b' = nullmat(`m_sd_1_b'), r(sd_1)
			mat `m_sd_2_b' = nullmat(`m_sd_2_b'), r(sd_2)
			mat `obs' = nullmat(`obs'), r(N_1)+r(N_2)
			mat `obs' = nullmat(`obs'), r(N_1)
			mat `obs' = nullmat(`obs'), r(N_2)
						
		qui reg `var' mother if round == 2010, cluster(id_mun)
		
			matrix b = e(b)
			matrix v = e(V)
			mat `m_coef_c_b' = nullmat(`m_coef_c_b'), b[1,1]
			mat `m_d_se_c_b' = nullmat(`m_d_se_c_b'), sqrt(v[1,1])
			mat `m_d_p_c_b' = nullmat(`m_d_p_c_b'), 2*ttail(e(df_r),abs(b[1,1]/sqrt(v[1,1])))
			mat `obs' = nullmat(`obs'), e(N)
					
		qui reg `var' mother $controls if round == 2010, cluster(id_mun)
		
			matrix b = e(b)
			matrix v = e(V)
			mat `cm_coef_c_b' = nullmat(`cm_coef_c_b'), b[1,1]
			mat `cm_d_se_c_b' = nullmat(`cm_d_se_c_b'), sqrt(v[1,1])
			mat `cm_d_p_c_b' = nullmat(`cm_d_p_c_b'), 2*ttail(e(df_r),abs(b[1,1]/sqrt(v[1,1])))
			mat `obs' = nullmat(`obs'), e(N)			
					
		qui ttest `var' if round == 2012, by(mother)
		
			mat `m_mu_f' = nullmat(`m_mu_f'), (r(mu_1)*r(N_1)+r(mu_2)*r(N_2))/(r(N_1)+r(N_2))
			mat `m_mu_1_f' = nullmat(`m_mu_1_f'), r(mu_1)
			mat `m_mu_2_f' = nullmat(`m_mu_2_f'), r(mu_2)
			mat `m_sd_f' = nullmat(`m_sd_f'), r(sd)
			mat `m_sd_1_f' = nullmat(`m_sd_1_f'), r(sd_1)
			mat `m_sd_2_f' = nullmat(`m_sd_2_f'), r(sd_2)
			mat `obs' = nullmat(`obs'), r(N_1)+r(N_2)
			mat `obs' = nullmat(`obs'), r(N_1)
			mat `obs' = nullmat(`obs'), r(N_2)
						
		qui reg `var' mother if round == 2012, cluster(id_mun)
		
			matrix b = e(b)
			matrix v = e(V)
			mat `m_coef_c_f' = nullmat(`m_coef_c_f'), b[1,1]
			mat `m_d_se_c_f' = nullmat(`m_d_se_c_f'), sqrt(v[1,1])
			mat `m_d_p_c_f' = nullmat(`m_d_p_c_f'), 2*ttail(e(df_r),abs(b[1,1]/sqrt(v[1,1])))
			mat `obs' = nullmat(`obs'), e(N)
			
		qui reg `var' mother $controls if round == 2012, cluster(id_mun)
		
			matrix b = e(b)
			matrix v = e(V)
			mat `cm_coef_c_f' = nullmat(`cm_coef_c_f'), b[1,1]
			mat `cm_d_se_c_f' = nullmat(`cm_d_se_c_f'), sqrt(v[1,1])
			mat `cm_d_p_c_f' = nullmat(`cm_d_p_c_f'), 2*ttail(e(df_r),abs(b[1,1]/sqrt(v[1,1])))
			mat `obs' = nullmat(`obs'), e(N)	
						
		qui treatreg `var' $controls if round == 2012, treat(mother_t = mother $controls) cluster(id_mun)
		
			matrix b = e(b)
			matrix v = e(V)
			mat `ivm_coef_c_f' = nullmat(`ivm_coef_c_f'), _b[mother_t]
			mat `ivm_d_se_c_f' = nullmat(`ivm_d_se_c_f'), _se[mother_t]
			mat `ivm_d_p_c_f' = nullmat(`ivm_d_p_c_f'), 2*ttail(e(N)-e(df_m)-1,abs(_b[mother_t]/_se[mother_t]))
			mat `obs' = nullmat(`obs'), e(N)	
			
	}
		
	foreach mat in m_mu_b m_mu_1_b m_mu_2_b m_sd_b m_sd_1_b m_sd_2_b m_coef_c_b m_d_se_c_b m_d_p_c_b cm_coef_c_b cm_d_se_c_b cm_d_p_c_b m_mu_f m_mu_1_f m_mu_2_f m_sd_f m_sd_1_f m_sd_2_f m_coef_c_f m_d_se_c_f m_d_p_c_f ivm_coef_c_f ivm_d_se_c_f ivm_d_p_c_f cm_coef_c_f cm_d_se_c_f cm_d_p_c_f {
		mat coln ``mat'' = `varlist'
	}

	tempname b V
	mat `b' = `m_coef_c_b'*0
	mat `V' = `b''*`b'
	eret post `b' `V'
	eret local cmd "treattest"
	
	mat temp = `obs'
	eret scalar N1 = temp[1,1] // Full sample baseline
	eret scalar N2 = temp[1,2] // HH sample baseline
	eret scalar N3 = temp[1,3] // Mother sample baseline
	eret scalar N4 = temp[1,4] // Regression sample baseline
	eret scalar N5 = temp[1,5] // Regression sample baseline with controls
	eret scalar N6 = temp[1,6] // Full sample follow-up
	eret scalar N7 = temp[1,7] // HH sample follow-up
	eret scalar N8 = temp[1,8] // Mother sample follow-up
	eret scalar N9 = temp[1,9] // Regression sample follow-up
	eret scalar N10 = temp[1,10] // Regression sample follow-up controls
	eret scalar N11 = temp[1,11] // IV Regression sample follow-up
		
	foreach mat in m_mu_b m_mu_1_b m_mu_2_b m_sd_b m_sd_1_b m_sd_2_b m_coef_c_b m_d_se_c_b m_d_p_c_b cm_coef_c_b cm_d_se_c_b cm_d_p_c_b m_mu_f m_mu_1_f m_mu_2_f m_sd_f m_sd_1_f m_sd_2_f m_coef_c_f m_d_se_c_f m_d_p_c_f ivm_coef_c_f ivm_d_se_c_f ivm_d_p_c_f cm_coef_c_f cm_d_se_c_f cm_d_p_c_f  obs {
		eret mat `mat' = ``mat''
	}
	
	end


********************************************************************************************
**** TREATTEST BOOTSTRAP
**** This program compares means and equality of means using bootstrap
********************************************************************************************

	capt prog drop treattest_bs

	program treattest_bs, eclass

	syntax varlist [if] [in], [ * ]
	marksample touse
	markout `touse'
	tempname m_mu_b m_mu_1_b m_mu_2_b m_sd_b m_sd_1_b m_sd_2_b m_coef_c_b m_d_se_c_b m_d_p_c_b m_coef_c_b_bs m_d_se_c_b_bs m_d_p_c_b_bs m_mu_f m_mu_1_f m_mu_2_f m_sd_f m_sd_1_f m_sd_2_f m_coef_c_f m_d_se_c_f m_d_p_c_f m_coef_c_f_bs m_d_se_c_f_bs m_d_p_c_f_bs obs                       
                           
	foreach var of local varlist {
		
		* Treatment 
				
		qui ttest `var' if `touse' & round == 2010, by(mother) `options'
		
			mat `m_mu_b' = nullmat(`m_mu_b'), (r(mu_1)*r(N_1)+r(mu_2)*r(N_2))/(r(N_1)+r(N_2))
			mat `m_mu_1_b' = nullmat(`m_mu_1_b'), r(mu_1)
			mat `m_mu_2_b' = nullmat(`m_mu_2_b'), r(mu_2)
			mat `m_sd_b' = nullmat(`m_sd_b'), r(sd)
			mat `m_sd_1_b' = nullmat(`m_sd_1_b'), r(sd_1)
			mat `m_sd_2_b' = nullmat(`m_sd_2_b'), r(sd_2)
			mat `obs' = nullmat(`obs'), r(N_1)+r(N_2)
			mat `obs' = nullmat(`obs'), r(N_1)
			mat `obs' = nullmat(`obs'), r(N_2)
						
		qui reg `var' mother if `touse' & round == 2010, cluster(id_mun)
		
			matrix b = e(b)
			matrix v = e(V)
			mat `m_coef_c_b' = nullmat(`m_coef_c_b'), b[1,1]
			mat `m_d_se_c_b' = nullmat(`m_d_se_c_b'), sqrt(v[1,1])
			mat `m_d_p_c_b' = nullmat(`m_d_p_c_b'), 2*ttail(e(df_r),abs(b[1,1]/sqrt(v[1,1])))
			mat `obs' = nullmat(`obs'), e(N)
		
		qui bootstrap, reps(500) cluster(id_mun): reg `var' mother if `touse' & round == 2010, cluster(id_mun)
		
			matrix b = e(b)
			matrix v = e(V)
			mat `m_coef_c_b_bs' = nullmat(`m_coef_c_b_bs'), b[1,1]
			mat `m_d_se_c_b_bs' = nullmat(`m_d_se_c_b_bs'), sqrt(v[1,1])
			mat `m_d_p_c_b_bs' = nullmat(`m_d_p_c_b_bs'), 2 * ttail(e(N)-e(df_m)-1, abs(b[1,1] / sqrt(v[1,1])))
			mat `obs' = nullmat(`obs'), e(N)
			
		qui ttest `var' if `touse' & round == 2012, by(mother) `options'
		
			mat `m_mu_f' = nullmat(`m_mu_f'), (r(mu_1)*r(N_1)+r(mu_2)*r(N_2))/(r(N_1)+r(N_2))
			mat `m_mu_1_f' = nullmat(`m_mu_1_f'), r(mu_1)
			mat `m_mu_2_f' = nullmat(`m_mu_2_f'), r(mu_2)
			mat `m_sd_f' = nullmat(`m_sd_f'), r(sd)
			mat `m_sd_1_f' = nullmat(`m_sd_1_f'), r(sd_1)
			mat `m_sd_2_f' = nullmat(`m_sd_2_f'), r(sd_2)
			mat `obs' = nullmat(`obs'), r(N_1)+r(N_2)
			mat `obs' = nullmat(`obs'), r(N_1)
			mat `obs' = nullmat(`obs'), r(N_2)
						
		qui reg `var' mother if `touse' & round == 2012, cluster(id_mun)
		
			matrix b = e(b)
			matrix v = e(V)
			mat `m_coef_c_f' = nullmat(`m_coef_c_f'), b[1,1]
			mat `m_d_se_c_f' = nullmat(`m_d_se_c_f'), sqrt(v[1,1])
			mat `m_d_p_c_f' = nullmat(`m_d_p_c_f'), 2*ttail(e(df_r),abs(b[1,1]/sqrt(v[1,1])))
			mat `obs' = nullmat(`obs'), e(N)			
		
		qui bootstrap, reps(500) cluster(id_mun): reg `var' mother if `touse' & round == 2012, cluster(id_mun)
		
			matrix b = e(b)
			matrix v = e(V)
			mat `m_coef_c_f_bs' = nullmat(`m_coef_c_f_bs'), b[1,1]
			mat `m_d_se_c_f_bs' = nullmat(`m_d_se_c_f_bs'), sqrt(v[1,1])
			mat `m_d_p_c_f_bs' = nullmat(`m_d_p_c_f_bs'), 2 * ttail(e(N)-e(df_m)-1, abs(b[1,1] / sqrt(v[1,1])))
			mat `obs' = nullmat(`obs'), e(N)	
						
	}
		
	foreach mat in m_mu_b m_mu_1_b m_mu_2_b m_sd_b m_sd_1_b m_sd_2_b m_coef_c_b m_d_se_c_b m_d_p_c_b m_coef_c_b_bs m_d_se_c_b_bs m_d_p_c_b_bs m_mu_f m_mu_1_f m_mu_2_f m_sd_f m_sd_1_f m_sd_2_f m_coef_c_f m_d_se_c_f m_d_p_c_f m_coef_c_f_bs m_d_se_c_f_bs m_d_p_c_f_bs {
		mat coln ``mat'' = `varlist'
	}

	tempname b V
	mat `b' = `m_coef_c_b'*0
	mat `V' = `b''*`b'
	eret post `b' `V'
	eret local cmd "treattest_bs"
	
	mat temp = `obs'
	eret scalar N1 = temp[1,1] // Full sample baseline
	eret scalar N2 = temp[1,2] // HH sample baseline
	eret scalar N3 = temp[1,3] // Mother sample baseline
	eret scalar N4 = temp[1,4] // Regression sample baseline
	eret scalar N5 = temp[1,5] // Regression sample baseline (bs)
	eret scalar N6 = temp[1,6] // Full sample follow-up
	eret scalar N7 = temp[1,7] // HH sample follow-up
	eret scalar N8 = temp[1,8] // Mother sample follow-up
	eret scalar N9 = temp[1,9] // Regression sample follow-up
	eret scalar N10 = temp[1,10] // Regression sample follow-up (bs)
		
	foreach mat in m_mu_b m_mu_1_b m_mu_2_b m_sd_b m_sd_1_b m_sd_2_b m_coef_c_b m_d_se_c_b m_d_p_c_b m_coef_c_b_bs m_d_se_c_b_bs m_d_p_c_b_bs m_mu_f m_mu_1_f m_mu_2_f m_sd_f m_sd_1_f m_sd_2_f m_coef_c_f m_d_se_c_f m_d_p_c_f m_coef_c_f_bs m_d_se_c_f_bs m_d_p_c_f_bs obs {
		eret mat `mat' = ``mat''
	}
	
	end


********************************************************************************************
**** TREATTEST CONTROLS (treattest_ct)
**** This program compares means and equality of means
********************************************************************************************

	capt prog drop treattest_ct

	program treattest_ct, eclass

	syntax varlist [if] [in], [ * ]
	marksample touse
	markout `touse'
	tempname m_coef_c_b1 m_d_se_c_b1 m_d_p_c_b1 m_coef_c_b2 m_d_se_c_b2 m_d_p_c_b2 m_coef_c_b3 m_d_se_c_b3 m_d_p_c_b3 m_coef_c_f1 m_d_se_c_f1 m_d_p_c_f1 m_coef_c_f2 m_d_se_c_f2 m_d_p_c_f2 m_coef_c_f3 m_d_se_c_f3 m_d_p_c_f3 obs
                                      
	foreach var of local varlist {
	
		bys id_mun round: egen `var'_base = mean(`var') if round == 2010

		bys id_mun: egen temp = max(`var'_base)
		replace `var'_base = temp
		drop temp
	
	
	}
	
	foreach var of local varlist {
		
		* BASELINE 
										
			qui reg `var' mother if `touse' & round == 2010, cluster(id_mun)
		
				matrix b = e(b)
				matrix v = e(V)
				mat `m_coef_c_b1' = nullmat(`m_coef_c_b1'), b[1,1]
				mat `m_d_se_c_b1' = nullmat(`m_d_se_c_b1'), sqrt(v[1,1])
				mat `m_d_p_c_b1' = nullmat(`m_d_p_c_b1'), 2 * ttail(e(df_r), abs(b[1,1] / sqrt(v[1,1])))
				mat `obs' = nullmat(`obs'), e(N)

			qui reg `var' mother $controls if `touse' & round == 2010, cluster(id_mun)
		
				matrix b = e(b)
				matrix v = e(V)
				mat `m_coef_c_b2' = nullmat(`m_coef_c_b2'), b[1,1]
				mat `m_d_se_c_b2' = nullmat(`m_d_se_c_b2'), sqrt(v[1,1])
				mat `m_d_p_c_b2' = nullmat(`m_d_p_c_b2'), 2 * ttail(e(df_r),abs(b[1,1] / sqrt(v[1,1])))
				mat `obs' = nullmat(`obs'), e(N)

			qui reg `var' mother $controls `var'_base if `touse' & round == 2010, cluster(id_mun)
		
				matrix b = e(b)
				matrix v = e(V)
				mat `m_coef_c_b3' = nullmat(`m_coef_c_b3'), b[1,1]
				mat `m_d_se_c_b3' = nullmat(`m_d_se_c_b3'), sqrt(v[1,1])
				mat `m_d_p_c_b3' = nullmat(`m_d_p_c_b3'), 2 * ttail(e(df_r),abs(b[1,1] / sqrt(v[1,1])))
				mat `obs' = nullmat(`obs'), e(N)
				
		* FOLLOW-UP
										
			qui reg `var' mother if `touse' & round == 2012, cluster(id_mun)
		
				matrix b = e(b)
				matrix v = e(V)
				mat `m_coef_c_f1' = nullmat(`m_coef_c_f1'), b[1,1]
				mat `m_d_se_c_f1' = nullmat(`m_d_se_c_f1'), sqrt(v[1,1])
				mat `m_d_p_c_f1' = nullmat(`m_d_p_c_f1'), 2 * ttail(e(df_r), abs(b[1,1] / sqrt(v[1,1])))
				mat `obs' = nullmat(`obs'), e(N)

			qui reg `var' mother $controls if `touse' & round == 2012, cluster(id_mun)
		
				matrix b = e(b)
				matrix v = e(V)
				mat `m_coef_c_f2' = nullmat(`m_coef_c_f2'), b[1,1]
				mat `m_d_se_c_f2' = nullmat(`m_d_se_c_f2'), sqrt(v[1,1])
				mat `m_d_p_c_f2' = nullmat(`m_d_p_c_f2'), 2 * ttail(e(df_r),abs(b[1,1] / sqrt(v[1,1])))
				mat `obs' = nullmat(`obs'), e(N)

			qui reg `var' mother $controls `var'_base if `touse' & round == 2012, cluster(id_mun)
		
				matrix b = e(b)
				matrix v = e(V)
				mat `m_coef_c_f3' = nullmat(`m_coef_c_f3'), b[1,1]
				mat `m_d_se_c_f3' = nullmat(`m_d_se_c_f3'), sqrt(v[1,1])
				mat `m_d_p_c_f3' = nullmat(`m_d_p_c_f3'), 2 * ttail(e(df_r),abs(b[1,1] / sqrt(v[1,1])))
				mat `obs' = nullmat(`obs'), e(N)
		
		* Drop average dependent variable at baseline
		
			drop `var'_base
								
	}
		
	foreach mat in m_coef_c_b1 m_d_se_c_b1 m_d_p_c_b1 m_coef_c_b2 m_d_se_c_b2 m_d_p_c_b2 m_coef_c_b3 m_d_se_c_b3 m_d_p_c_b3 m_coef_c_f1 m_d_se_c_f1 m_d_p_c_f1 m_coef_c_f2 m_d_se_c_f2 m_d_p_c_f2 m_coef_c_f3 m_d_se_c_f3 m_d_p_c_f3 {
		mat coln ``mat'' = `varlist'
	}

	tempname b V
	mat `b' = `m_coef_c_b1'*0
	mat `V' = `b''*`b'
	eret post `b' `V'
	eret local cmd "treattest_ct"
	
	mat temp = `obs'
	eret scalar N1 = temp[1,1] 
	eret scalar N2 = temp[1,2] 
	eret scalar N3 = temp[1,3] 
	eret scalar N4 = temp[1,4] 
	eret scalar N5 = temp[1,5] 
	eret scalar N6 = temp[1,6] 
		
	foreach mat in m_coef_c_b1 m_d_se_c_b1 m_d_p_c_b1 m_coef_c_b2 m_d_se_c_b2 m_d_p_c_b2 m_coef_c_b3 m_d_se_c_b3 m_d_p_c_b3 m_coef_c_f1 m_d_se_c_f1 m_d_p_c_f1 m_coef_c_f2 m_d_se_c_f2 m_d_p_c_f2 m_coef_c_f3 m_d_se_c_f3 m_d_p_c_f3 obs {
		eret mat `mat' = ``mat''
	}
	
	end



********************************************************************************************
**** SUM_EXP
**** This program decomposes expenditure and consumption for different items
********************************************************************************************

capt prog drop sum_exp

	program sum_exp, eclass

	syntax anything [if] [in], [ * ]
	marksample touse
	markout `touse'
	tempname m_mu1 m_sd1 m_mu2 m_sd2 m_mu3 m_sd3 m_mu4 m_sd4 m_mu5 m_sd5 m_mu6 m_sd6 obs
                           
	foreach x of local anything {
		
	* Expenditure shares
				
		qui ttest sh_`x' if `touse', by(mother) `options'
		
			mat `m_mu1' = nullmat(`m_mu1'), (r(mu_1)*r(N_1)+r(mu_2)*r(N_2))/(r(N_1)+r(N_2))
			mat `m_sd1' = nullmat(`m_sd1'), r(sd)
			mat `obs' = nullmat(`obs'), r(N_1)+r(N_2)
		
	* Consumption shares
		
		qui ttest sh_`x'_adj if `touse', by(mother) `options'
		
			mat `m_mu2' = nullmat(`m_mu2'), (r(mu_1)*r(N_1)+r(mu_2)*r(N_2))/(r(N_1)+r(N_2))
			mat `m_sd2' = nullmat(`m_sd2'), r(sd)
			mat `obs' = nullmat(`obs'), r(N_1)+r(N_2)

	* Decomposition of consumption shares
	
		qui ttest sh_pur_`x'_adj if `touse', by(mother) `options'
		
			mat `m_mu3' = nullmat(`m_mu3'), (r(mu_1)*r(N_1)+r(mu_2)*r(N_2))/(r(N_1)+r(N_2))
			mat `m_sd3' = nullmat(`m_sd3'), r(sd)
			mat `obs' = nullmat(`obs'), r(N_1)+r(N_2)

		qui ttest sh_sf_`x'_adj if `touse', by(mother) `options'
		
			mat `m_mu4' = nullmat(`m_mu4'), (r(mu_1)*r(N_1)+r(mu_2)*r(N_2))/(r(N_1)+r(N_2))
			mat `m_sd4' = nullmat(`m_sd4'), r(sd)
			mat `obs' = nullmat(`obs'), r(N_1)+r(N_2)

		qui ttest sh_stor_`x'_adj if `touse', by(mother) `options'
		
			mat `m_mu5' = nullmat(`m_mu5'), (r(mu_1)*r(N_1)+r(mu_2)*r(N_2))/(r(N_1)+r(N_2))
			mat `m_sd5' = nullmat(`m_sd5'), r(sd)
			mat `obs' = nullmat(`obs'), r(N_1)+r(N_2)

		qui ttest sh_gift_`x'_adj if `touse', by(mother) `options'
		
			mat `m_mu6' = nullmat(`m_mu6'), (r(mu_1)*r(N_1)+r(mu_2)*r(N_2))/(r(N_1)+r(N_2))
			mat `m_sd6' = nullmat(`m_sd6'), r(sd)
			mat `obs' = nullmat(`obs'), r(N_1)+r(N_2)


		
	}
	
	foreach mat in m_mu1 m_sd1 m_mu2 m_sd2 m_mu3 m_sd3 m_mu4 m_sd4 m_mu5 m_sd5 m_mu6 m_sd6 obs {
		mat coln ``mat'' = `anything'
	}
	
	tempname b V
	mat `b' = `m_mu1'*0
	mat `V' = `b''*`b'
	eret post `b' `V'
	eret local cmd "sum_exp"
	
	mat temp = `obs'
	eret scalar N1 = temp[1,1]
	eret scalar N2 = temp[1,2]
	eret scalar N3 = temp[1,3]
	eret scalar N4 = temp[1,4]
	eret scalar N5 = temp[1,5]
	eret scalar N6 = temp[1,6]
		
	foreach mat in m_mu1 m_sd1 m_mu2 m_sd2 m_mu3 m_sd3 m_mu4 m_sd4 m_mu5 m_sd5 m_mu6 m_sd6 obs {
		eret mat `mat' = ``mat''
	
	}
	
	end	

